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We study entanglement creation between two independent XX chains, which are repeatedly cou¬ 
pled locally to spin-1/2 Bell pairs. We show analytically that in the steady state the entanglement 
of the Bell pairs is perfectly transferred to the chains, generating large-scale interchain pair cor¬ 
relations. However, before the steady state is reached, within a growing causal region around the 
interacting locus the chains are found in a current driven nonquilibrium steady state (NESS). In 
the NESS, the chains cross entanglement decays exponentially with respect to the distance to the 
boundary sites with a typical length scale which is inversely proportional to the driving current. 

PACS numbers: 03.67.Bg, 03.65.Yz, 42.50.Dv 


Entanglement plays one of the major roles on the scene 
of quantum engineering and quantum control [Tj 2]. This 
purely quantum property is critical in quantum informa¬ 
tion processing [3H3- Hence entanglement is the center 
of attention of numerous studies ISHD1- In a recent work, 
S. Zippilli and co-workers studied the steady state entan¬ 
glement replication in two independent quantum many- 
body systems, both locally driven by a common entan¬ 
gled field 0 flTj. They showed that the field can be tuned 
to perfectly replicate the driving entanglement across the 
initially independent arrays through the generation of a 
scale-free set of two-particle Bell states. 

In this work, we study the dynamics of entanglement 
replication between two A'X-quantum spin chains each 
of size N. The chains are both locally driven at one 
of their boundaries by an entangled Markovian environ¬ 
ment. Starting from any Gaussian initial state, the sys¬ 
tem reaches, for t tn (with tat ~ N 3 ), a unique sta¬ 
tionary state (already observed in [8] [PL]) in which spins 
belonging to different chains are perfectly entangled by 
pairs. During the relaxation toward the stationary state, 
the system passes through a transient regime for inter¬ 
mediate times Td <C t •C N/v c , where is a typical mi¬ 
croscopic relaxation time related to the local dissipative 
coupling to the bath and v c the typical sound velocity 
within the chains. In this time regime, close to the inter¬ 
acting boundaries the chains fall into a Non-Equilibrium 
Steady State (NESS) with a steady current j z * injected 
by the environment into the chains. In the NESS, the 
entanglement between spins belonging to different chains 
and facing each other decay exponentially over a length 
scale £ as we move away from the interacting boundaries. 
This entanglement length scale £ is shown to be propor¬ 
tional to the inverse of the steady current j z *. Moreover, 
as the NESS extends ballistically through the chains, lo¬ 
cal observables along the chains show typically a scaling 
form Q(x,t) — fc}{x/t) where x is the distance from the 
interacting boundaries. 

The system is composed of two non-interacting XX 



FIG. 1: (Color online) Two non interacting XX chains are 
coupled on one edge, via a RIP, to one component of a Bell 
pair. The first chain is labeled from 1 to N, and the second 
one from N + 3 to 2N + 2 

chains with Hamiltonian 

H s = -^~ X («+i + 0-K+1) - | X W 

n6 S* ne s 

with S = {1,2,... ,N} U {iV + 3, N + 4,..., 2N + 2} 
and S* = S\{N,2N + 2}. The hopping constant K is 
fixed to 1/2 such that the sound velocity v c = 1 into 
the chains. Note that the total transverse magnetization 
M z = Yjn a n commutes with Hs . As a consequence, 
we will show that the external held does not affect the 
dynamics of the system, when considering Gaussian ini¬ 
tial state. The Markovian environment with which the 
system interacts is described by a Repeated Interaction 
Process (RIP) [I2HI6] . see Fig. 0 for a pictorial repre¬ 
sentation. More precisely, the environment is made of an 
infinite set of identical and independent spin pairs, with 
Hamiltonians h[^ = —h{al x + crjr; 2 )/2 where k € IN*, 

such that the bath Hamiltonian is Hb = Y^kLi^B^■ 
Each bath spin pair is prepared into a perfectly entan¬ 
gled Bell state r? = |$)($| with |$) = ^ (| U) + | 4-t», 
where | f) and I i) are eigenstates of the erjF Pauli matri- 
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ces. Therefore, the initial state of the bath is given by the 
direct product t)b = 0 fcgN „ V- During the RIP the Bell 
pairs r)k will interact one after the other, over a typical 
interaction time t, with the XX-chains. The dynamics of 
the full system + environment is generated by the time 
dependent Hamiltonian H(t) = Hs + Hb + V(t), where 
V ( t ) describes the system-environment coupling and is 
given by V(t) = V for t £ [(fc — l)r, kr[. The inter¬ 
action couples one spin of the pair, say to the left 
chain only while the second one, 0 ^ 2 , is only coupled 
to the right one and that through a local isotropic XX 
coupling 

y(k) = ~ 9 (°^ CT M + a N aV k,l + a h^N +3 + ^,2^+3) • 

( 2 ) 

Note that the coupling amplitude 7 is taken to be the 
same for the left and right chains. 

The initial system + environment state is taken to be 
the tensor product state p( 0 ) = qs(0)®Vb where the sys¬ 
tem density matrix gg( 0 ) is a fully factorized (thermal) 
mixture £>s( 0 ) = (£) neS g n with the single spin density 
matrix g n = (t n + g, n a^)/2, where fj, n = Tr{a^g n }. At 
g, n = 1 (— 1 ) the density matrix g n reduces to the pure 
state | t)(t 1(1 l)U !)• 

Over one time step of the RIP, the two chains, to¬ 
gether with the Bell pair with which they interact, for¬ 
mally constitute a XX chain of size 2 N + 2 with non- 
homogeneous couplings. This chain can be mapped to 
a free Fermi system through a Jordan Wigner transfor¬ 
mation c n = YlU(-^)°n [E] with 2 cr± = a*± icr y . 
The c and adjoint c' satisfy the canonical anticommu¬ 
tation relations { c\,Cj} = 5ij and {cj,cj} = {cj,Cj} = 
0. To proceed further we separate the system vari¬ 
ables from the bath ones. We define the field operator 

'P t = (xt.ft) with xt = (4 ,...,c] v ,c] v+ 3,...,4jv+ 2 ) 

and ft = (cjy , 1; c' N+2 )- The total Hamiltonian takes the 
form H = where T is a (2N + 2) x (21V+ 2) matrix 

given by 


T = 


(T S 7©\ 

V701 TbJ 


— ht, 


T S 



(3) 


where A (of size N x N) is given by Aij = —K(6ij +1 + 
Sij-i) and 0 (of size 2 N x 2 ) is given by 0, :j = 
— (D.tvU, 1 + As long as we are not interested 

in the fate of the environment Bell pair, it is unneces¬ 
sary to further specify the 2x2 matrix Tb as it does not 
appear in the dynamical equation which described the 
system evolution [14] . 

The initial state of the ensemble system+environment 
being Gaussian [22] . the total density matrix g{t) and 
the system density matrix gs(t) = Trslgit)} remain 
Gaussian during the time evolution generated by the 
free fermionic Hamiltonian [18] . Thanks to Wick the¬ 
orem, a full description of the system is given in terms 
of its two points correlators only. As a first step we de¬ 
fine the ( 2N + 2) x (2 N + 2) correlations matrix G by 


Gij(t) = (T,|\l !j)(t). Its evolution over one time step is 
given by G(r) = e~ lrT G{0)e lrT . The magnetic field h 
appearing on the diagonal of T only, the previous equa¬ 
tion clearly indicates that h has no effect on the system 
dynamics. Next we define the reduced 2N x 2N correla¬ 
tions matrix Gs by (Gs)ij(t) = (x]xj)(t). When taking 
the continuous time limit (r —>• 0 ), one has to keep q 2 r 
constant i 12]. This procedure leads to a differential equa¬ 
tion describing the evolution of the system in interaction 
with the environment. The system correlation matrix 
obeys a Lindblad-like differential equation El 

d t Gs = -C(G S ) + r 2 0G B 0 t , (4) 

with £(.) = i[Ts, .]+T 2 {., 00f}/2 where T is the rescaled 
coupling parameter. The 2x2 matrix Gb is the corre¬ 
lation matrix of a bath pair, explicitly given in the Bell 
state ?7 by ( Gs)ij = 1/2 Vi, j. Choosing the coupling pa¬ 
rameter T = 7 y/r leads to the approximate solution for 
the discrete case. Corrections of order 1/r are, however 
expected when considering the discrete case. 

In the following we use the concurrence as a measure 
of the entanglement between two spins k, l £ S of the 
system. For an arbitrary density matrix g 1 representing 
the quantum state of the two spins, the concurrence m 
1201 is defined by C = 2max{0, Ai — A 2 — A 3 — A 4 } where the 
A’s are the square roots of the eigenvalues (in descending 
order) of the matrix R = gg with g = (a y <£>a v )g(cr v ®<j y ) 
and where g is the complex conjugate of g (in the a z 
diagonal base). 

We start, with the simplest case where the chains are 
reduced to two qubits. In such a case Tg = 0, 0 = — 12 x 2 
and Eq. 0 becomes d t Gs = —T 2 (Gs — Gb), with solu¬ 
tion (c{c 4 )(f) = (l-e -r “f)/ 2 , (CjCj)(t) = (l+g,je~ r ~t)/2 
for j = 1,4. From there, we can reconstruct the reduced 
density matrix 4 1 ’ 4 ) and evaluate the concurrence: 

C^(t) = max jo, 1 - e- r2 ‘ — \ x (5) 

\f [2e~ r2t + (mfu - l)e“ 2r24 ]" - (/ui + iM) 2 e~ 2r2t j . 

Interestingly enough, one observes a finite delay time T, 
defined as the waiting time needed to generate entangle¬ 
ment [tV 1 , 4 )(f) = 0 for all t <T\, for initially mixed states 
(for fi n ^ ±1). On the contrary, when the two spins are 
initially prepared in a pure state then T vanishes. As 
an example, in the case of opposite initial magnetiza¬ 
tions /ii = — ^ = g one has T = T -2 In 

which vanishes as (1 — l^l ) 1 ^ 2 f° r —i► ±1. This fi¬ 
nite waiting time T observed for mixed initial states 
can be interpreted as the time needed for quantum cor¬ 
relations to overcome the thermal fluctuations, before 
the entanglement can start to grow more or less with 
the same functional dependence as in the pure initial 
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state case. At long times, independently of the sys¬ 
tem initial state the concurrence C( 1,4 \t —► oo) = 1 
reflecting the fact that the system correlation matrix 
Gg(t —► oo) = Gg. Consequently, in the steady state, 
the non vanishing system’s density matrix coefficients are 
given by gy = gy = - 1 / 2 , gy = gy = 1/2 cor¬ 
responding to the Bell state |$~) = ^(| fl) — | It))- In 
other words, the RIP leads to a transfer of the entangle¬ 
ment from the reservoir to the system. 

We consider now the most interesting case of two 
chains of equal size N. The chains are supposed to be pre¬ 
pared into the same factorized homogeneous state, with 
initial magnetization g on each site. The hopping con¬ 
stant K in (JTJ) is fixed to 1 /2 such that the sound velocity 
v, : = 1 into the chains. 

Transient NESS: For sufficiently large system sizes N 
and times t < N, close to the boundaries directly in 
contact with the environment the system behaves like 
a semi-infinite one: no excitations have enough time to 
travel from one boundary to the other. In this regime, 
the system approaches locally a NESS with a finite flux 
of excitations injected at the contact point with the en¬ 
vironment and traveling into the chains away from the 
system-environment interaction site. This flux of par¬ 
ticles propagates the new state into the chains, modi¬ 
fying the values of the local quantities and eventually 
generating correlations. During this process, the trans- 
versed magnetization converges algebraically towards a 
constant value m z *: m z (p,t ) — m z * oc t~ 2 , as illustrated 
in Fig. [2]- (a). Outside the immediate vicinity of the 
boundary point p = 0 (which corresponds to the sites 
directly in contact with the bath), the deviation from 
the asymptotic flat stationary magnetization profile takes 
the scaling form m z (p,t) — m z * = f m {p/t ), as clearly 
seen on Fig. §( c). The local magnetization current 
j z (p,t)i which is defined through the continuity equa¬ 
tion d t m z (x,t) + d x j z (x,t) = 0 decays as t~ 3 toward an 
asymptotic value j z * as shown in Fig. ®b) for different 
locus p along the chains. Since the current j z (p,t) sat¬ 
isfy the continuity equation, it takes also a scaling form 
j z (p,t) = fj(p/t) as seen on Fig. |2[(d). 

For t < N, the repeated interactions with the envi¬ 
ronment generate in the chains steady longitudinal and 
cross entanglement in the vicinity of the contact point 
within a region which grows linearly in time. This en¬ 
tanglement is measured respectively by the longitudinal 
Ci(p) and cross C c (p) concurrences. The later is defined 
between two spins facing each other (one on site N — p 
and the other on N + 3 + p as parametrized in Fig. [lj. 
The former, is measured between two neighboring sites 
p and p + 1 on a given chain. Numerical results ob¬ 
tained from exact diagonalization (see | 2 T] for details) 
are shown in Fig. [3] for g = 1. Mimicking the behavior 
of the local (in chain) quantities like the magnetization 
and the current, the longitudinal concurrence Ci [p, t) con- 



FIG. 2: (Color online) Time evolution of the magnetization 
and current for different positions for two chains of size N = 
300, T = 0.5 with initial magnetization equal to 1. 

verges toward a steady value by following a scaling form 
Ci(p,t) = fciip/t )> which is clearly visible on Fig. |3[(b). 
Along the repeated interaction process, as already stated 
cross entanglement is generated between the two chains 
measured by C c (p,t) and which converges algebraically 
toward a steady value C c (p). For g = 1 the steady cross 
concurrence C c {p) decays exponentially with the distance 
p from the contact point as seen in Figs. [3]-a) and[3]-(c): 

Gc(p) oc exp (-p/£) . ( 6 ) 

The decay is over a typical length scale £(T) which seems 
to be inversely proportional to the steady current j 2 *(r) 
as seen in Fig. [3]-(d). The relation £(T) a l/j 2 *(F) is 
particularly good for small and large values of T. Devia¬ 
tions can be observed for values around F = 1. Note the 
symmetry j z *( F) = j z *(l/T) as in [14] . 

For 0 < g < 1 [23], the behavior of C c (p) is a bit more 
involved and it is reported in Fig. [4] for two different 
values of T, different positions p and as a function of 
g. We see that the cross concurrence decreases, until it 
reaches 0, as the initial magnetization g decreases. The 
value of g bellow which a given pair p disentangles con¬ 
verges rapidly toward the saturation value g = 1 as the 
pair label p is increased. This means that at finite initial 
temperature (g < 1 ) very few pairs are entangled in the 
NESS. For T = 1/2 we also see that there is a threshold 
value gt.hre bellow which there is no more cross entan¬ 
glement. On the contrary at larger dissipation rate F, 
for example T = 2 , the first pair is always entangled to a 
very high value (close to one) whatever the initial magne¬ 
tization is, the price to be paid being that the remaining 
pairs are already disentangled for values of g < 0.9. 

This intriguing behavior of the cross concurrence can 
be explained in the following way (at least for the first 
pair p = 0 ): assuming that the magnetization of spins 
N and N + 3 in the NESS is proportional to the initial 
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FIG. 3: (Color online) (a): Concurrence in the pairs p as a 
function of time. The size of the chains is N = 60 and the 
system reservoir coupling is T = 0.5. (b) Longitudinal con¬ 
currence profile at different times as a function of the rescaled 
parameter p/t of a chain of size N = 500 and T = 0.5. (c) 
Logarithm of the pair concurrence at t = N as a function of 
the pair p for three different values of T and N = 60. (d) 
Inverse of the stationary current j z * and entanglement length 
£ taken at t = N as a function of the logarithm of the dissi¬ 
pation coupling T for two chains of size N = 60. /3 ~ 5.67 is 
a proportionality coefficient. For all plots, the initial magne¬ 
tization of the chains is p = 1. 



FIG. 4: (Color online) (a) Cross concurrence in different pairs 
as a function of the initial magnetization p for F = 1/2. (b) 
Same as (a) for T = 2. The size of the chains is N = 60. 
The dots are the numerical results and the dashed line is the 
theoretical prediction. 

magnetization of the chain, itin,n +3 = (1 — oc)p and that 
the fermionic correlator is proportional to the Bell pairs 
correlations, (c^Cat+ 3) = a(cj v+1 c/v+2 ) = ck/ 2, where 
the coefficient a is a function of T. It follows by solv¬ 
ing locally the steady Lindblad equations that the cross 
concurrence of the first pair p = 0 is given by 

C c (0) = max jo, a - /r)| (7) 

with g(a,p) = (1 — a) 2 (l — /rt 2 ) [(1 + a) 2 — p 2 (l — a) 2 ]. 
For p = 1 one has C c ( 0) = a which fixes then the value 
of a for a given T. Depending on the value of a there 
may or may not exist a threshold value of p, pthre defined 
by C c ( 0) = 0 \/p < pthrei below which the entanglement 


in the first pair is lost in the NESS. For a > \[2 — 1 
(which corresponds to T ~ 0.5916), the first pair is always 
entangled in the NESS whatever the value of p. In Fig. 
[4] the theoretical prediction 0 is plotted for the two 
different values of T and one can see that it matches 
perfectly the numerical results. 

After the NESS regime, the system undergoes a re¬ 
laxation regime which results from the superposition of 
constructive and destructive interferences of excitations 
traveling up and down the chains, being reflected by the 
boundaries many times. This process leads to an ex¬ 
ponential relaxation towards the stationary state with a 
relaxation time diverging with the system size (r oc N 3 ). 
In the following we focus on the stationary state reached 
in the limit t —>■ 00 and satisfying C(G* S ) = r 2 0Gs0 1 '. 

Steady state. At very large times, due to multiple 
reflections at the chains boundaries of the excitations 
injected by the bath, the local NESS is progressively 
lost. The system converges toward its real unique sta¬ 
tionary state with a typical relaxation time of order N 3 . 
This final stationary state can be obtained in the fol¬ 
lowing way: Split the correlation function Gs,b into a 
real and imaginary part, Gk = Nk + iJk (k = S, B) 
with Nj£ = Ni- and = — J^. The J5 matrix con¬ 
tains current-like terms of the form i(c' n c m — c' m c n ) which 
have to vanish in the steady state and the remaining 
matrix Ns has to satisfy from Q the steady equations 
[T S ,N S ] = 0 , {00t,A/c?} = 20iV B 0t. The only non¬ 
vanishing entries of the solution matrix Ns are the on¬ 
site densities (cjc,-) = 1/2 and the cross-correlations 
(c\c 2 N+ 3 -i) = 1/2, Vi € S. It follows that only sites fac¬ 
ing each other are correlated and the steady state turns 
out to be given by an alternation of pairs in the per¬ 
fectly entangled states l^) = -^=(| fl) ± | It))- The 
steady density matrix factorizes into g* s = (iSipQp with 
Q 2 n = |$~>($ - | and Q 2 n +1 = |-F+)(4>+|. 

In summary the replication mechanism is perfect in the 
steady state. However, and most interestingly, in an in¬ 
termediate NESS regime, for pure initial state (\p\ = 1), 
the cross entanglement between the chains decays expo¬ 
nentially with respect to the distance to the interacting 
locus. This exponential decay is over a typical length 
scale which is found to be inversely proportional to the 
driving current. At finite temperatures (\p\ < 1) and at 
low enough dissipation rate, there is a threshold magne¬ 
tization below which the cross entanglement between the 
chains is completely lost in the NESS. 
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